#%% [markdown]
""" 
现给到北京线上某款视频APP在2018年1月-2018年7月的实际用户数据，需要预测2018年8月-2019年12月的数据，请基于课程中的案例分析，建立分析方案并给到预测的数据结果。
"""
# %%
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy import optimize
from statsmodels.formula import api as sfm

plt.rcParams["font.family"] = "Microsoft YaHei"
#%%
# 数据导入和转换
ldf = pd.read_csv("task2.csv", encoding="gbk")
df = pd.DataFrame(data={"temp": [i for i in range(24)]})
df = pd.merge(
    df,
    ldf,
    how="left",
    left_index=True,
    right_index=True,
)
df["year"] = 2018 + df["temp"] // 12
df["month"] = df["temp"] % 12 + 1
df["date"] = df["year"].astype(str) + "-" + df["month"].astype(str)
df.drop(columns=["temp", "日期"], inplace=True)
df.rename(columns={"当月实际户数（户）": "users"}, inplace=True)
#%%
# 产品上线后可能存在活动或峰值,需要考虑如何排除
# 通过对怎张率筛选,如果增长率超过其几何平均值,则认为当前处于活动中
df["cha"] = df.users.diff(1)
ldf["ch_rate"] = (df.users - df.users.shift(1)) / df.users.shift(1)
df["ch_rate"] = (df.users - df.users.shift(1)) / df.users.shift(1)
df["is_event"] = (df.users - df.users.shift(1)) / df.users.shift(1) > np.product(
    ldf.ch_rate
) ** (1 / len(ldf.ch_rate))
df
#%%
# 图表绘制
plt.figure(figsize=(12, 5))
plt.scatter(x=df.index, y=df.users)
plt.xticks(df.index, df.date, rotation=-90)
plt.yticks([i * 500 for i in range(7)])


# %% [markdown]
""" 
模型选择
+ 自回归
+ 线性模型
+ 非线性模型
    + 指数
    + 多项式
    + 对数
评估指标
"""
#%%
# 评估指标
def _statistics(y, y_prodict):
    t_res = stats.ttest_rel(y, y_prodict)
    f_res = stats.f_oneway(y, y_prodict)
    ssr = np.sum((y_prodict - np.mean(y)) ** 2)
    sse = np.sum((y - y_prodict) ** 2)
    sst = np.sum((y - np.mean(y)) ** 2)
    r2 = 1 - sse / sst

    MSE = np.mean(np.square(y - y_prodict))
    RMSE = np.sqrt(np.mean(np.square(y - y_prodict)))
    MAE = np.mean(np.abs(y - y_prodict))
    MAPE = np.mean(np.abs((y - y_prodict) / y)) * 100
    return {
        "t_test": t_res,
        "ANOVA": f_res,
        "r^2": {"r^2": r2, "sse": sse, "sst": sst, "ssr": ssr},
        "other": {"MSE": MSE, "RMSE": RMSE, "MAE": MAE, "MAPE": MAPE},
    }


#%%
def linerfit(model, y_prodict, s=""):
    plt.figure(figsize=(12, 5))
    plt.scatter(x=df.index, y=df.users)
    plt.plot(t.index, y_prodict, c="r", linestyle=":", marker="o")
    plt.xticks(df.index, df.date, rotation=-90)
    ticker_r = [i * 500 for i in range(7)]
    plt.yticks([i * 500 for i in range(7)])
    plt.text(0, max(ticker_r) / 2, s=s.format(model), fontsize=14)
    return _statistics(t.users, y_prodict)


t = df[["users", "year", "is_event"]].dropna()
t.reset_index(inplace=True)
#%% [markdown]
""" 
三元哑变量回归未能通过模型的假设检验,多个指标t检验未能通过
进行向后剔除
"""
#%%
model = sfm.ols(formula="users~index^3+index^2+index+is_event", data=t).fit()
model.summary()
#%% [markdown]
""" 
是否参与活动任然未能通过
进行向后剔除
"""
#%%
model = sfm.ols(formula="users~index^3+index+is_event", data=t).fit()
model.summary()
#%% [markdown]
""" 
各项指标接近通过,只有截距未能通过
绘图检验后发现模型偏离非常严重,不建议使用
"""
#%%
model = sfm.ols(formula="users~index^3+index", data=t).fit()
model.summary()
#%%
t["prodict"] = (
    model.params[0] + model.params[1] * t.index ** 3 + model.params[2] * t.index
)
linerfit(
    model.params, t.prodict, s="users={0[0]: .2f} + {0[1]: .2f}index^3+{0[2]: .2f}index"
)
plt.figure(figsize=(12, 5))
plt.scatter(x=t.index, y=t.users - t["prodict"])
#%%
model = sfm.ols(formula="users~index^2+index", data=t).fit()
model.summary()
#%%
t["prodict"] = (
    model.params[0] + model.params[1] * t.index ** 2 + model.params[2] * t.index
)
plt.figure(figsize=(12, 5))
plt.scatter(x=t.index, y=t.users - t["prodict"])
linerfit(model.params, t["prodict"], "一元y={0[0]:.2f}+{0[1]:.2f}x^2+{0[2]:.2f}x")
#%%
model = sfm.ols(formula="users~index/is_event", data=t).fit()
model.summary()
#%%
model = sfm.ols(formula="users~index + is_event", data=t).fit()
model.summary()
#%%
model = sfm.ols(formula="users~index", data=t).fit()
model.summary()
t["prodict"] = model.params[0] + model.params[1] * t.index
plt.figure(figsize=(12, 5))
plt.scatter(x=t.index, y=t.users - t["prodict"])
linerfit(model.params, t["prodict"], "一元y={0[0]:.2f}+{0[1]:.2f}x")

# %%
def exp_func(x, a, b):
    return a * np.exp(-b * x)


exp_model, modle_cov = optimize.curve_fit(exp_func, t.index, t.users)
exp_model
exp_model, modle_cov
y_prodict = np.array([exp_func(i, exp_model[0], exp_model[1]) for i in range(7)])
linerfit(
    exp_model,
    y_prodict,
    s="指数 y={0[0]:.2f}e^({0[1]:.2f}x)",
)

plt.figure(figsize=(12, 5))
plt.scatter(x=t.index, y=t.users - y_prodict)

# %% [markdown]
""" 
+ 经过测试
    + 发现还是一般线性回归拟合效果较好
    + 虽然多项式回归(3次)的拟合度更高,但是由于属于强行拟合因此后续会造成较大偏差,不与使用
    + 怎张率并不稳定所以也无法作为参考列
    + 哑变量回归 is_event 未能通过模型检验,且R^2提升较少,因此也已于剔除模型
    + 残差分布情况,线性回归的残差分布相较于其他模型更为集中,效果更好
+ 最终确认模型采用
+ $ user_prodict = 115.7143 + 133.4286x $

{'t_test': Ttest_relResult(statistic=-8.031425095730656e-16, pvalue=0.9999999999999993),
 'ANOVA': F_onewayResult(statistic=2.797963436614649e-31, pvalue=1.0),
 'r^2': {'r^2': 0.9388450031210315,
  'sse': 32470.85714285715,
  'sst': 530960.0,
  'ssr': 498489.14285714296},
 'other': {'MSE': 4638.693877551022,
  'RMSE': 68.10795752003595,
  'MAE': 63.42857142857144,
  'MAPE': 16.804138516363036}}

"""

# %%
model = sfm.ols(formula="users~index", data=t).fit()
model.summary()
#%%
df["prodict"] = model.params[0] + model.params[1] * df.index
df
# %%
plt.figure(figsize=(12, 5))
plt.plot(df.index, df.prodict, label="预测值", c="r")
plt.text(
    0,
    df.prodict.max() / 2,
    s="一元y={0[0]:.2f}x+{0[1]:.2f}".format(model.params),
    fontsize=14,
    c="r",
)
plt.scatter(df.index, df.users)
plt.figure(figsize=(12, 5))
plt.scatter(x=t.index, y=t.users - t["prodict"])
# %% [markdown]
""" 
最终值
"""
#%%
df[['date','users','prodict']]
# %%
